library(knitr)
options(max.print="75")
opts_chunk$set(echo = TRUE, cache = TRUE, warning = FALSE, prompt=FALSE, tidy = TRUE, comment = NA, message = FALSE)
opts_knit$set(width=75)

# Loading library here
library(tidyverse)
library(ggpubr)
library(MASS)
library(segmented)
library(agricolae)
library(car)
library(effectsize)
library(rstatix)
library(patchwork)

# set path
source("assets/utils.R")

Initial setup

Setup color index

# Exposure level
Expo_color <- c("#1b9e77", "#d95f02", "#8da0cb")
names(Expo_color) <- c("Low", "Medium", "High")
print(Expo_color)
      Low    Medium      High 
"#1b9e77" "#d95f02" "#8da0cb" 
#plot(1:3, 1:3, col = Expo_color, pch = 16, cex = 4)
# Ethnicity
Ethnicity_color <- gg_color_hue(3)
names(Ethnicity_color) <- c("SANEMA", "YEKWANA", "Visitors")
print(Ethnicity_color)
   SANEMA   YEKWANA  Visitors 
"#F8766D" "#00BA38" "#619CFF" 
Sanemas <- c("Chajuranha", "Mosenahanha", "Kuyuwininha", "Shianana-Jiyakwanha", "Washudihanha", "Sudukuma")
Yekwana <- c("Kanarakuni", "Fiyakwanha")

Import the data

# Alpha diversities
alpha = readRDS("../data/alpha_imported_all.RDS")
# Mapping
## This mapping files include all villagers and baseline of visitors of all body sites
mt_ms = readRDS("../data/Mapping_MS_Expo.Rds")
Mapping_work <- mt_ms %>% mutate(Expo = ifelse(Ethnicity=="SANEMA", "Low", ifelse(Village=="Fiyakwanha", "Low", ifelse(SampleGroup=="Visitors", "High", "Medium"))) %>% factor(levels = c("Low", "Medium", "High")), Age_num = as.numeric(Age), Age_num = ifelse(is.na(Age_num), as.numeric(gsub(pattern = "_months", replacement = "", Age, ignore.case = T))/12, Age_num), Age_grp1 = ifelse(Age_num>=18, "Adults", "Children"))
Mapping_work$Age_num[Mapping_work$Age=="4_Days"] <- 0
Mapping_work$Age_num[Mapping_work$Age=="1_day"] <- 0

Set up the Analysis

dat_alpha <- merge(alpha, Mapping_work, by = 1) %>% mutate(Age_grp2 = case_when(Age_num<=3 ~ "Age_0-3", Age_num<=8 ~ "Age_3-8", Age_num<18 ~ "Age_8-18", T ~ "Adults"))

Effects of exposure level

Indices <- colnames(alpha[c(2:5)])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Mouth", "Nose", "Right_Arm", "Right_Hand")

Initial test on feces (no need to run)

We first check the faith_pd which described the overall complexity of the microbial community.

BB = "Feces"
AA = "faith_pd"
# Year 2015
YY = "2015"
Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==BB, Year==YY, SampleGroup=="Villagers")

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender)

anova(regforward.out)

The best model from 2015 includes Expo, Age, Expo:Age interaction and Gender (in decreasing importance).

# Year 2016
YY = "2016"
Work_dat <- Work %>% dplyr::filter(Body_Site==BB, Year==YY, SampleGroup=="Villagers")

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender)

anova(regforward.out)

The best model from 2016 includes Expo, Age, Expo:Age interaction (in decreasing importance), but Gender was not there.

So we will fitting the data with Expo, Age, Expo:Age interaction, and Gender for both years to keep unity.

# Final model
lm_models <- list()
for (YY in Yrs){
  Work_dat <- Work %>% dplyr::filter(Body_Site==BB, Year==YY, SampleGroup=="Villagers")
  fit.final <- lm(as.formula(paste0(AA, " ~ Age_num + Expo + Age_num:Expo + Gender")), data = Work_dat)
  print(paste0("--------Below is ANOVA table for year ", YY))
  anova(fit.final) %>% print()
  
  print(paste0("--------Below is coefficient for year ", YY))
  summary(fit.final) %>% print()
  lm_models[[YY]] <- fit.final
}

Next we perform the broken stick regression based on the selected final model

for (YY in Yrs){
  #YY = "2016"
  fit.seg <- segmented(lm_models[[YY]], seg.Z = ~ Age_num, psi = c(3))
  summary(fit.seg)
  plot(fit.seg$fitted.values, fit.seg$residuals)
  plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
  plot(fit.seg$model$Age_num, hatvalues(fit.seg))
  plot(fit.seg$model$Age_num, cooks.distance(fit.seg))
  plot(hatvalues(fit.seg), fit.seg$residuals)
}

It seems from the residual plots that the point with faith_pd smaller than 10 is dragging the curve down and causing non-random residual.

Work %>% dplyr::filter(Body_Site==BB, faith_pd<12, SampleGroup=="Villagers")

It does appears that many baby are outliers, so re-do the analysis without baby.

Modeling without babies

Overview

for (bb in BodySites){
    for (aa in Indices){
        Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==bb, SampleGroup=="Villagers", Age_num>=1)
        p <- ggplot(Work_dat, aes(x = Age_num, y = get(aa), color = Expo, shape = Gender)) +
            geom_point(size = 1.5) +
            scale_shape_manual(values = c(16, 21)) +
            scale_color_manual(breaks = names(Expo_color), values = Expo_color) +
            facet_grid(.~Year) +
            theme_jw() + theme(aspect.ratio = 0.8, panel.background = element_rect(color = "black")) +
            labs(x = "Age (Years)", y = aa, title = bb)
        print(p)
    }
}

#### Exploratory Statisticals ##### Children Building full multiple linear regression models with interaction terms for children

alpha_youngchild_lm_models = list()
lm.stats.tbl = list()
for (bb in BodySites){
    for (aa in Indices){
        for (yy in Yrs){
            Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==bb, Year==yy, SampleGroup=="Villagers", Age_num>=1, Age_num<18)
            lm.model <- lm(paste0(aa, " ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender"), data = Work_dat, contrasts = list(Expo = "contr.sum", Gender = "contr.sum"))
            lm.stats = Anova(lm.model, type = 3) %>% # Type 3 ss
                merge(., effectsize::eta_squared(., ci = 0.95), by.x = 0, by.y = 1, all = T) %>% # effect size estimates
                rename(Factors = Row.names) %>% mutate(Year = yy, Body_Site = bb, Indice = aa)
            alpha_youngchild_lm_models[[bb]][[aa]][[yy]] = lm.model
            lm.stats.tbl[[paste(bb, aa, yy, sep = "-")]] = lm.stats
        }
    }
}
lm.stats.tbls = do.call("rbind", lm.stats.tbl) %>% rename(p_value = `Pr(>F)`) %>% mutate(q_value = p.adjust(p_value, method = "fdr"))

# full model (multiple linear regression) statistics for alpha diversity in children
#write.table(lm.stats.tbls, file = "../output/tables/alpha_lm_stats_children.txt", sep = "\t", row.names = F, quote = F)

lm.stats.tbls.2 = lm.stats.tbls %>% dplyr::filter(!is.na(Eta2_partial)) %>% mutate(Factors = factor(Factors, levels = c("Age_num", "Gender", "Expo", "Age_num:Expo", "Age_num:Gender", "Expo:Gender"))) %>% arrange(Body_Site, Indice, Year, Factors)

# p values    
ggplot(lm.stats.tbls.2, aes(x = Factors, y = q_value, color = Body_Site)) +
    geom_point() +
    geom_path(aes(group = Year, linetype = Year)) +
    facet_grid(Body_Site~Indice) +
    scale_y_continuous(breaks = c(0, 0.05, 0.1, 0.2, 1.0)) +
    theme_jw(base_size = 12) + theme(aspect.ratio = 0.4, axis.text.x = element_text(angle = 90, hjust = 1), panel.grid.minor = element_blank()) +
    ggtitle("FDR adjusted p-value for full models of alpha diversity in Children")

# eta partial
ggplot(lm.stats.tbls.2, aes(x = Factors, y = Eta2_partial, color = Body_Site)) +
    geom_point() +
    geom_path(aes(group = Year, linetype = Year)) +
    facet_grid(Body_Site~Indice, scales = "free_y") +
    #scale_y_continuous(breaks = c(0, 0.05, 0.1, 0.2, 1.0)) +
    theme_jw() + theme(aspect.ratio = 0.4, axis.text.x = element_text(angle = 90, hjust = 1)) +
    ggtitle("Effect size (eta partial) for full models of alpha diversity in Children")

No strong interactions obsereved among children

all ages

Building full multiple linear regression models with interaction terms

alpha_full_lm_models = list()
lm.stats.tbl = list()
for (bb in BodySites){
    for (aa in Indices){
        for (yy in Yrs){
            Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==bb, Year==yy, SampleGroup=="Villagers", Age_num>=1)
            lm.model <- lm(paste0(aa, " ~ Age_num + Expo + Gender + Age_num:Expo + Age_num:Gender + Expo:Gender"), data = Work_dat, contrasts = list(Expo = "contr.sum", Gender = "contr.sum"))
            lm.stats = Anova(lm.model, type = 3) %>% # Type 3 ss
                merge(., effectsize::eta_squared(., ci = 0.95), by.x = 0, by.y = 1, all = T) %>% # effect size estimates
                rename(Factors = Row.names) %>% mutate(Year = yy, Body_Site = bb, Indice = aa)
            alpha_full_lm_models[[bb]][[aa]][[yy]] = lm.model
            lm.stats.tbl[[paste(bb, aa, yy, sep = "-")]] = lm.stats
        }
    }
}
lm.stats.tbl = do.call("rbind", lm.stats.tbl)
#o full model (multiple linear regression) statistics for alpha diversity in all age
write.table(lm.stats.tbl, file = "../results/tables/alpha_lm_stats.txt", sep = "\t", row.names = F, quote = F)

lm.stats.tbl.2 = lm.stats.tbl %>% dplyr::filter(!is.na(Eta2_partial)) %>% rename(p_value = `Pr(>F)`) %>% mutate(Factors = factor(Factors, levels = c("Age_num", "Gender", "Expo", "Age_num:Expo", "Age_num:Gender", "Expo:Gender"))) %>% arrange(Body_Site, Indice, Year, Factors)
    
ggplot(lm.stats.tbl.2, aes(x = Factors, y = p_value, color = Body_Site)) +
    geom_point() +
    geom_path(aes(group = Year, linetype = Year)) +
    facet_grid(Body_Site~Indice) +
    scale_y_continuous(breaks = c(0, 0.05, 0.1, 0.2, 1.0)) +
    coord_cartesian(ylim = c(0, 0.2)) +
    theme_jw() + theme(aspect.ratio = 0.6, axis.text.x = element_text(angle = 90, hjust = 1)) +
    guides(color = F) +
    ggtitle("FDR adjusted p-value for full models of alpha diversity in all age")

ggplot(lm.stats.tbl.2, aes(x = Factors, y = Eta2_partial, color = Body_Site)) +
    geom_point() +
    geom_path(aes(group = Year, linetype = Year)) +
    facet_grid(Body_Site~Indice, scales = "free_y") +
    theme_jw() + theme(aspect.ratio = 0.6, axis.text.x = element_text(angle = 90, hjust = 1)) +
    guides(color = F) +
    ggtitle("Effect size (eta-partial) for full models of alpha diversity in all age")

#p1 + p2 + plot_layout(guides = "collect")

Final Models (No interaction terms)

models = list()
model_tbls = data.frame(BodySite = as.character(),
                        Index = as.character(),
                        Year = as.character(),
                        Model = as.character(),
                        Model.type = as.character())

for (bb in BodySites){
    for (aa in Indices){
        for (yy in Yrs){
            Work_dat <- dat_alpha %>% dplyr::filter(Body_Site==bb, Year==yy, SampleGroup=="Villagers", Age_num>=1)
            # Obtain the full linear model
            fit.full <- lm(paste0(aa, " ~ Age_num + Expo + Gender"), data = Work_dat)
            # Obtain the null model with no predictors, this is necessary for forward selection use stepAIC
            fit.null <- lm(paste0(aa, " ~ 1"), data = Work_dat)
            
            # Forward selection, scope parameters provides the full model
            regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~ Age_num + Expo + Gender)
            
            lm.model <- lm(regforward.out$terms, data = Work_dat, contrasts = list(Expo = "contr.sum", Gender = "contr.sum"))
            lm.model.anova = car::Anova(lm.model, type = 3)
            
            if ("Age_num" %in% rownames(lm.model.anova)){
                fit.seg <- segmented(regforward.out, seg.Z = ~ Age_num, psi = 9)
                model_comp = anova(fit.seg, lm.model)
                if (!is.na(model_comp[2, 6]) & model_comp[2, 6]<=0.05){
                    fit.final = fit.seg
                    summary(fit.seg) %>% print()
                    davies.test(lm.model, seg.Z = ~ Age_num, k = 10) %>% print()
                    confint.segmented(fit.seg) %>% print()
                    slope(fit.seg) %>% print()
                } else {
                    fit.final = lm.model
                }
            } else {
                fit.final = lm.model
            }
            plot(fit.final$fitted.values, fit.final$residuals)
            plot((fit.final$model)[, 1], fit.final$fitted.values)
            plot(hatvalues(fit.final), fit.final$residuals)
            qqnorm(fit.final$residuals)
            ff = (fit.final$terms %>% as.character())[3]
            models[[bb]][[aa]][[yy]] = fit.final
            if ("segmented" %in% class(fit.final)){
                model_tbls = rbind(model_tbls, cbind(bb, aa, yy, ff, "segmented"))
            } else {
                model_tbls = rbind(model_tbls, cbind(bb, aa, yy, ff, "linear"))
            }
        }
    }
}
Start:  AIC=577.4
faith_pd ~ 1

          Df Sum of Sq    RSS    AIC
+ Expo     1    877.60 5581.9 556.91
+ Age_num  1    251.81 6207.6 573.27
+ Gender   1    126.68 6332.8 576.35
<none>                 6459.5 577.40

Step:  AIC=556.91
faith_pd ~ Expo

          Df Sum of Sq    RSS    AIC
+ Gender   1    177.43 5404.4 553.94
+ Age_num  1    134.92 5446.9 555.14
<none>                 5581.9 556.91

Step:  AIC=553.94
faith_pd ~ Expo + Gender

          Df Sum of Sq    RSS    AIC
+ Age_num  1    164.54 5239.9 551.17
<none>                 5404.4 553.94

Step:  AIC=551.17
faith_pd ~ Expo + Gender + Age_num


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
               Est. St.Err
psi1.Age_num 7.918   1.69

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  20.2457     2.0376   9.936  < 2e-16 ***
ExpoMedium   -4.3660     0.9399  -4.645 7.45e-06 ***
GenderM      -2.2691     0.9107  -2.492  0.01382 *  
Age_num       1.1480     0.4147   2.768  0.00636 ** 
U1.Age_num   -1.1658     0.4157  -2.804       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.578 on 148 degrees of freedom
Multiple R-Squared: 0.287,  Adjusted R-squared: 0.2629 

Convergence attained in 3 iter. (rel. change 1.9747e-16)

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Expo + Gender + Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.0001387
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 7.91839     4.57902    11.2578
$Age_num
            Est.  St.Err.  t value CI(95%).l CI(95%).u
slope1  1.148000 0.414700  2.76820  0.328470  1.967500
slope2 -0.017832 0.035963 -0.49584 -0.088899  0.053236

Start:  AIC=610.01
faith_pd ~ 1

          Df Sum of Sq    RSS    AIC
+ Expo     1    466.10 5026.9 596.32
+ Age_num  1     74.17 5418.9 609.60
<none>                 5493.0 610.01
+ Gender   1     34.72 5458.3 610.89

Step:  AIC=596.32
faith_pd ~ Expo

          Df Sum of Sq    RSS    AIC
<none>                 5026.9 596.32
+ Gender   1    55.767 4971.2 596.34
+ Age_num  1    44.091 4982.8 596.76

Start:  AIC=-798.63
pielou_e ~ 1

          Df Sum of Sq     RSS     AIC
<none>                 0.85047 -798.63
+ Expo     1 0.0072315 0.84324 -797.95
+ Gender   1 0.0067885 0.84368 -797.87
+ Age_num  1 0.0013146 0.84915 -796.87

Start:  AIC=-928.53
pielou_e ~ 1

          Df Sum of Sq     RSS     AIC
+ Expo     1  0.048001 0.87411 -935.99
<none>                 0.92211 -928.53
+ Age_num  1  0.000622 0.92149 -926.65
+ Gender   1  0.000009 0.92210 -926.53

Step:  AIC=-935.99
pielou_e ~ Expo

          Df  Sum of Sq     RSS     AIC
<none>                  0.87411 -935.99
+ Gender   1 0.00035701 0.87375 -934.07
+ Age_num  1 0.00002161 0.87409 -934.00

Start:  AIC=-52.93
shannon ~ 1

          Df Sum of Sq    RSS     AIC
+ Expo     1    5.8680 101.93 -59.555
+ Gender   1    1.9180 105.88 -53.699
<none>                 107.80 -52.935
+ Age_num  1    0.7776 107.02 -52.050

Step:  AIC=-59.55
shannon ~ Expo

          Df Sum of Sq    RSS     AIC
+ Gender   1   2.41770  99.51 -61.251
<none>                 101.93 -59.555
+ Age_num  1   0.28031 101.65 -57.979

Step:  AIC=-61.25
shannon ~ Expo + Gender

          Df Sum of Sq    RSS     AIC
<none>                 99.510 -61.251
+ Age_num  1   0.44589 99.064 -59.943

Start:  AIC=-99.94
shannon ~ 1

          Df Sum of Sq    RSS      AIC
+ Expo     1    7.3360 92.171 -111.493
<none>                 99.507  -99.938
+ Age_num  1    0.3042 99.203  -98.480
+ Gender   1    0.0459 99.461  -98.019

Step:  AIC=-111.49
shannon ~ Expo

          Df Sum of Sq    RSS     AIC
<none>                 92.171 -111.49
+ Gender   1  0.168447 92.003 -109.82
+ Age_num  1  0.091109 92.080 -109.67

Start:  AIC=1366.98
observed_otus ~ 1

          Df Sum of Sq     RSS    AIC
+ Expo     1    109046  979633 1352.7
+ Age_num  1     41607 1047072 1363.0
+ Gender   1     21417 1067263 1365.9
<none>                 1088679 1367.0

Step:  AIC=1352.73
observed_otus ~ Expo

          Df Sum of Sq    RSS    AIC
+ Gender   1     28701 950932 1350.2
+ Age_num  1     24608 955025 1350.8
<none>                 979633 1352.7

Step:  AIC=1350.15
observed_otus ~ Expo + Gender

          Df Sum of Sq    RSS    AIC
+ Age_num  1     29694 921238 1347.3
<none>                 950932 1350.2

Step:  AIC=1347.26
observed_otus ~ Expo + Gender + Age_num


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
               Est. St.Err
psi1.Age_num 8.333  1.747

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  170.512     24.597   6.932 1.19e-10 ***
ExpoMedium   -48.451     12.635  -3.835 0.000186 ***
GenderM      -29.088     12.246  -2.375 0.018818 *  
Age_num       12.481      4.284   2.913 0.004131 ** 
U1.Age_num   -12.605      4.307  -2.927       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 75.27 on 148 degrees of freedom
Multiple R-Squared: 0.2298,  Adjusted R-squared: 0.2038 

Convergence attained in 4 iter. (rel. change 0)

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Expo + Gender + Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.001613
alternative hypothesis: two.sided

               Est. CI(95%).low CI(95%).up
psi1.Age_num 8.3333     4.88077    11.7858
$Age_num
           Est. St.Err.  t value CI(95%).l CI(95%).u
slope1 12.48100 4.28410  2.91330    4.0150  20.94700
slope2 -0.12465 0.52149 -0.23902   -1.1552   0.90588

Start:  AIC=1452.83
observed_otus ~ 1

          Df Sum of Sq    RSS    AIC
+ Expo     1     34581 607797 1445.0
<none>                 642379 1452.8
+ Age_num  1      3124 639254 1454.0
+ Gender   1      1626 640753 1454.4

Step:  AIC=1445.04
observed_otus ~ Expo

          Df Sum of Sq    RSS    AIC
<none>                 607797 1445.0
+ Gender   1    2901.4 604896 1446.2
+ Age_num  1    1507.6 606290 1446.6

Start:  AIC=503.39
faith_pd ~ 1

          Df Sum of Sq    RSS    AIC
+ Expo     1    327.20 3160.9 489.23
+ Age_num  1    326.20 3161.9 489.29
<none>                 3488.1 503.39
+ Gender   1      3.18 3484.9 505.24

Step:  AIC=489.23
faith_pd ~ Expo

          Df Sum of Sq    RSS    AIC
+ Age_num  1   218.027 2942.9 479.51
<none>                 3160.9 489.23
+ Gender   1     0.781 3160.1 491.19

Step:  AIC=479.51
faith_pd ~ Expo + Age_num

         Df Sum of Sq    RSS    AIC
<none>                2942.9 479.51
+ Gender  1  0.062452 2942.8 481.51

Start:  AIC=448.86
faith_pd ~ 1

          Df Sum of Sq    RSS    AIC
+ Age_num  1    578.15 2494.4 420.21
+ Expo     1    224.27 2848.3 439.71
<none>                 3072.6 448.86
+ Gender   1     31.84 3040.7 449.32

Step:  AIC=420.21
faith_pd ~ Age_num

         Df Sum of Sq    RSS    AIC
+ Expo    1   172.536 2321.9 411.68
<none>                2494.4 420.21
+ Gender  1     3.114 2491.3 422.03

Step:  AIC=411.68
faith_pd ~ Age_num + Expo

         Df Sum of Sq    RSS    AIC
<none>                2321.9 411.68
+ Gender  1   0.87557 2321.0 413.62

Start:  AIC=-736.73
pielou_e ~ 1

          Df Sum of Sq    RSS     AIC
+ Expo     1   0.33382 1.4799 -768.09
<none>                 1.8138 -736.73
+ Gender   1   0.01466 1.7991 -736.06
+ Age_num  1   0.00116 1.8126 -734.84

Step:  AIC=-768.09
pielou_e ~ Expo

          Df Sum of Sq    RSS     AIC
+ Gender   1 0.0225098 1.4574 -768.60
<none>                 1.4799 -768.09
+ Age_num  1 0.0067515 1.4732 -766.84

Step:  AIC=-768.6
pielou_e ~ Expo + Gender

          Df Sum of Sq    RSS     AIC
<none>                 1.4574 -768.60
+ Age_num  1 0.0050221 1.4524 -767.17

Start:  AIC=-674.18
pielou_e ~ 1

          Df Sum of Sq    RSS     AIC
<none>                 1.4779 -674.18
+ Expo     1 0.0086101 1.4692 -673.04
+ Gender   1 0.0045496 1.4733 -672.63
+ Age_num  1 0.0041333 1.4737 -672.59

Start:  AIC=-16.6
shannon ~ 1

          Df Sum of Sq    RSS     AIC
+ Expo     1   30.7103 115.70 -53.211
<none>                 146.41 -16.604
+ Age_num  1    1.1004 145.31 -15.841
+ Gender   1    0.3191 146.09 -14.962

Step:  AIC=-53.21
shannon ~ Expo

          Df Sum of Sq    RSS     AIC
<none>                 115.70 -53.211
+ Gender   1   0.70883 114.99 -52.219
+ Age_num  1   0.00265 115.70 -51.215

Start:  AIC=-10.31
shannon ~ 1

          Df Sum of Sq    RSS      AIC
+ Age_num  1    3.6509 131.54 -12.3362
<none>                 135.19 -10.3118
+ Gender   1    0.4367 134.75  -8.7875
+ Expo     1    0.2287 134.96  -8.5607

Step:  AIC=-12.34
shannon ~ Age_num

         Df Sum of Sq    RSS     AIC
<none>                131.54 -12.336
+ Gender  1   0.12639 131.41 -10.478
+ Expo    1   0.10877 131.43 -10.458

Start:  AIC=1202.36
observed_otus ~ 1

          Df Sum of Sq    RSS    AIC
+ Expo     1   25991.7 221496 1186.2
+ Age_num  1    9449.1 238038 1198.0
<none>                 247488 1202.4
+ Gender   1     191.2 247296 1204.2

Step:  AIC=1186.16
observed_otus ~ Expo

          Df Sum of Sq    RSS    AIC
+ Age_num  1    4428.0 217068 1184.8
<none>                 221496 1186.2
+ Gender   1      33.7 221462 1188.1

Step:  AIC=1184.85
observed_otus ~ Expo + Age_num

         Df Sum of Sq    RSS    AIC
<none>                217068 1184.8
+ Gender  1   0.49699 217067 1186.8

Start:  AIC=1099.97
observed_otus ~ 1

          Df Sum of Sq    RSS    AIC
+ Age_num  1   28638.6 229086 1084.7
+ Expo     1   13366.0 244359 1094.2
<none>                 257725 1100.0
+ Gender   1     416.3 257309 1101.7

Step:  AIC=1084.66
observed_otus ~ Age_num

         Df Sum of Sq    RSS    AIC
+ Expo    1   10545.4 218541 1079.7
<none>                229086 1084.7
+ Gender  1      51.1 229035 1086.6

Step:  AIC=1079.73
observed_otus ~ Age_num + Expo

         Df Sum of Sq    RSS    AIC
<none>                218541 1079.7
+ Gender  1    186.96 218354 1081.6

Start:  AIC=717.39
faith_pd ~ 1

          Df Sum of Sq   RSS    AIC
+ Expo     1   1204.32 24178 712.73
+ Age_num  1    784.55 24598 715.09
<none>                 25383 717.39
+ Gender   1      1.37 25381 719.39

Step:  AIC=712.73
faith_pd ~ Expo

          Df Sum of Sq   RSS    AIC
+ Age_num  1   1193.82 22985 707.80
<none>                 24178 712.73
+ Gender   1      0.87 24178 714.73

Step:  AIC=707.8
faith_pd ~ Expo + Age_num

         Df Sum of Sq   RSS    AIC
<none>                22985 707.80
+ Gender  1    13.127 22972 709.72

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 10.503  1.915

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.2593     3.5611   0.354 0.724183    
ExpoMedium    7.4982     2.1762   3.446 0.000764 ***
Age_num       2.0584     0.5730   3.592 0.000461 ***
U1.Age_num   -2.1586     0.5797  -3.724       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.12 on 132 degrees of freedom
Multiple R-Squared: 0.2359,  Adjusted R-squared: 0.2127 

Convergence attained in 3 iter. (rel. change 0)

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Expo + Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 7.96e-05
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 10.5028     6.71477    14.2907
$Age_num
           Est.  St.Err. t value CI(95%).l CI(95%).u
slope1  2.05840 0.573020  3.5922   0.92492  3.191900
slope2 -0.10021 0.096118 -1.0426  -0.29034  0.089921

Start:  AIC=661.41
faith_pd ~ 1

          Df Sum of Sq   RSS    AIC
+ Age_num  1    797.06 15612 656.54
+ Expo     1    287.61 16122 660.97
<none>                 16409 661.41
+ Gender   1     14.96 16394 663.29

Step:  AIC=656.54
faith_pd ~ Age_num

         Df Sum of Sq   RSS    AIC
+ Expo    1    515.66 15097 653.91
<none>                15612 656.54
+ Gender  1     82.15 15530 657.81

Step:  AIC=653.91
faith_pd ~ Age_num + Expo

         Df Sum of Sq   RSS    AIC
<none>                15097 653.91
+ Gender  1    64.915 15032 655.31

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   13  2.365

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.3916     2.8886   1.867 0.064172 .  
Age_num       1.3496     0.3750   3.599 0.000449 ***
ExpoMedium    4.6226     1.8022   2.565 0.011427 *  
U1.Age_num   -1.4254     0.3823  -3.729       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.854 on 133 degrees of freedom
Multiple R-Squared: 0.213,  Adjusted R-squared: 0.1894 

Convergence attained in 5 iter. (rel. change 3.6307e-07)

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Age_num + Expo ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0002563
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 12.9999     8.32219    17.6777
$Age_num
            Est.  St.Err.  t value CI(95%).l CI(95%).u
slope1  1.349600 0.374950  3.59950   0.60799  2.091300
slope2 -0.075799 0.082967 -0.91361  -0.23990  0.088305

Start:  AIC=-563.47
pielou_e ~ 1

          Df Sum of Sq    RSS     AIC
+ Expo     1  0.040152 2.1687 -563.98
<none>                 2.2088 -563.47
+ Age_num  1  0.017469 2.1914 -562.56
+ Gender   1  0.003786 2.2050 -561.70

Step:  AIC=-563.98
pielou_e ~ Expo

          Df Sum of Sq    RSS     AIC
<none>                 2.1687 -563.98
+ Age_num  1 0.0097853 2.1589 -562.60
+ Gender   1 0.0036188 2.1651 -562.21

Start:  AIC=-586.78
pielou_e ~ 1

          Df Sum of Sq    RSS     AIC
+ Expo     1  0.043216 1.8930 -587.89
<none>                 1.9363 -586.78
+ Gender   1  0.008933 1.9273 -585.41
+ Age_num  1  0.006969 1.9293 -585.27

Step:  AIC=-587.89
pielou_e ~ Expo

          Df Sum of Sq    RSS     AIC
<none>                 1.8930 -587.89
+ Gender   1 0.0122879 1.8808 -586.79
+ Age_num  1 0.0020227 1.8910 -586.04

Start:  AIC=67.68
shannon ~ 1

          Df Sum of Sq    RSS    AIC
+ Age_num  1   11.5234 209.75 62.353
<none>                 221.27 67.680
+ Gender   1    0.2747 221.00 69.510
+ Expo     1    0.2379 221.03 69.533

Step:  AIC=62.35
shannon ~ Age_num

         Df Sum of Sq    RSS    AIC
<none>                209.75 62.353
+ Gender  1   0.93295 208.82 63.742
+ Expo    1   0.01050 209.74 64.346

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 12.022  2.763

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.94666    0.31238   6.232 5.66e-09 ***
Age_num      0.14656    0.05043   2.906  0.00429 ** 
U1.Age_num  -0.15464    0.05135  -3.012       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.182 on 133 degrees of freedom
Multiple R-Squared: 0.1598,  Adjusted R-squared: 0.1409 

Convergence attained in 3 iter. (rel. change 3.0031e-07)

    Davies' test for a change in the slope

data:  formula = shannon ~ Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 0.0006814
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 12.0223     6.55758     17.487
$Age_num
             Est.   St.Err.  t value CI(95%).l CI(95%).u
slope1  0.1465600 0.0504330  2.90610  0.046809  0.246320
slope2 -0.0080792 0.0096185 -0.83996 -0.027104  0.010946

Start:  AIC=40.36
shannon ~ 1

          Df Sum of Sq    RSS    AIC
+ Age_num  1    7.5644 174.65 36.506
<none>                 182.22 40.357
+ Expo     1    0.4060 181.81 42.049
+ Gender   1    0.1408 182.08 42.250

Step:  AIC=36.51
shannon ~ Age_num

         Df Sum of Sq    RSS    AIC
<none>                174.65 36.506
+ Gender  1   0.77672 173.88 37.891
+ Expo    1   0.01420 174.64 38.494

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 17.227  3.354

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.17025    0.24915   8.710 9.96e-15 ***
Age_num      0.09580    0.02856   3.355  0.00103 ** 
U1.Age_num  -0.11278    0.03039  -3.711       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.069 on 134 degrees of freedom
Multiple R-Squared: 0.1591,  Adjusted R-squared: 0.1403 

Convergence attained in 4 iter. (rel. change 0)

    Davies' test for a change in the slope

data:  formula = shannon ~ Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0002751
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 17.2273     10.5927    23.8619
$Age_num
            Est.  St.Err. t value CI(95%).l CI(95%).u
slope1  0.095802 0.028556  3.3549  0.039324 0.1522800
slope2 -0.016981 0.010403 -1.6322 -0.037557 0.0035952

Start:  AIC=1264.61
observed_otus ~ 1

          Df Sum of Sq     RSS    AIC
+ Expo     1     45099 1332819 1262.0
<none>                 1377918 1264.6
+ Age_num  1     19166 1358752 1264.7
+ Gender   1       758 1377160 1266.5

Step:  AIC=1262.05
observed_otus ~ Expo

          Df Sum of Sq     RSS    AIC
+ Age_num  1     31672 1301147 1260.8
<none>                 1332819 1262.0
+ Gender   1       840 1331979 1264.0

Step:  AIC=1260.75
observed_otus ~ Expo + Age_num

         Df Sum of Sq     RSS    AIC
<none>                1301147 1260.8
+ Gender  1    2782.7 1298365 1262.5

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   10   1.86

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -14.592     27.044  -0.540 0.590405    
ExpoMedium    45.538     16.526   2.755 0.006690 ** 
Age_num       15.130      4.352   3.477 0.000687 ***
U1.Age_num   -16.118      4.402  -3.661       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 92.05 on 132 degrees of freedom
Multiple R-Squared: 0.1882,  Adjusted R-squared: 0.1636 

Convergence attained in 6 iter. (rel. change 1.8623e-07)

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Expo + Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.0002517
alternative hypothesis: two.sided

             Est. CI(95%).low CI(95%).up
psi1.Age_num   10      6.3204    13.6796
$Age_num
         Est. St.Err. t value CI(95%).l CI(95%).u
slope1 15.130 4.35170  3.4769    6.5225   23.7380
slope2 -0.988 0.72994 -1.3535   -2.4319    0.4559

Start:  AIC=1188.48
observed_otus ~ 1

          Df Sum of Sq    RSS    AIC
+ Age_num  1   24477.4 723352 1185.9
<none>                 747829 1188.5
+ Expo     1    4861.7 742967 1189.6
+ Gender   1    1190.4 746639 1190.3

Step:  AIC=1185.89
observed_otus ~ Age_num

         Df Sum of Sq    RSS    AIC
<none>                723352 1185.9
+ Expo    1     10229 713122 1185.9
+ Gender  1      4033 719319 1187.1

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   13  2.513

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   28.101     17.880   1.572 0.118390    
Age_num        8.146      2.338   3.485 0.000666 ***
U1.Age_num    -8.876      2.416  -3.675       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 69.02 on 134 degrees of freedom
Multiple R-Squared: 0.1463,  Adjusted R-squared: 0.1272 

Convergence attained in 6 iter. (rel. change 7.0875e-06)

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0009512
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 13.0001     8.03041    17.9698
$Age_num
           Est. St.Err. t value CI(95%).l CI(95%).u
slope1  8.14620 2.33770  3.4848    3.5227  12.77000
slope2 -0.73013 0.60844 -1.2000   -1.9335   0.47326

Start:  AIC=898.45
faith_pd ~ 1

          Df Sum of Sq   RSS    AIC
+ Expo     1    9453.0 29340 854.64
+ Gender   1     750.3 38043 897.24
<none>                 38793 898.45
+ Age_num  1     301.5 38492 899.17

Step:  AIC=854.64
faith_pd ~ Expo

          Df Sum of Sq   RSS    AIC
+ Gender   1    978.22 28362 851.08
<none>                 29340 854.64
+ Age_num  1      0.12 29340 856.64

Step:  AIC=851.08
faith_pd ~ Expo + Gender

          Df Sum of Sq   RSS    AIC
<none>                 28362 851.08
+ Age_num  1    1.1997 28361 853.07

Start:  AIC=457.7
faith_pd ~ 1

          Df Sum of Sq   RSS    AIC
+ Expo     1   1468.41 18646 453.41
<none>                 20115 457.70
+ Gender   1     61.22 20054 459.45
+ Age_num  1     51.33 20064 459.49

Step:  AIC=453.41
faith_pd ~ Expo

          Df Sum of Sq   RSS    AIC
<none>                 18646 453.41
+ Age_num  1   107.687 18539 454.93
+ Gender   1    52.232 18594 455.18

Start:  AIC=-905.64
pielou_e ~ 1

          Df Sum of Sq     RSS     AIC
+ Gender   1 0.0292817 0.61829 -911.23
<none>                 0.64757 -905.64
+ Expo     1 0.0052636 0.64231 -904.98
+ Age_num  1 0.0004897 0.64708 -903.76

Step:  AIC=-911.23
pielou_e ~ Gender

          Df Sum of Sq     RSS     AIC
<none>                 0.61829 -911.23
+ Expo     1 0.0063057 0.61199 -910.91
+ Age_num  1 0.0009688 0.61732 -909.49

Start:  AIC=-453.82
pielou_e ~ 1

          Df Sum of Sq     RSS     AIC
+ Expo     1  0.040233 0.30175 -462.21
+ Age_num  1  0.011223 0.33076 -454.59
<none>                 0.34199 -453.82
+ Gender   1  0.001064 0.34092 -452.08

Step:  AIC=-462.21
pielou_e ~ Expo

          Df Sum of Sq     RSS     AIC
+ Age_num  1 0.0151271 0.28663 -464.48
<none>                 0.30175 -462.21
+ Gender   1 0.0012784 0.30047 -460.56

Step:  AIC=-464.48
pielou_e ~ Expo + Age_num

         Df Sum of Sq     RSS     AIC
<none>                0.28663 -464.48
+ Gender  1 0.0016638 0.28496 -462.96

Start:  AIC=-67.4
shannon ~ 1

          Df Sum of Sq    RSS     AIC
+ Expo     1    6.6969 100.72 -75.954
+ Gender   1    4.8961 102.52 -73.048
<none>                 107.42 -67.397
+ Age_num  1    0.2545 107.16 -65.786

Step:  AIC=-75.95
shannon ~ Expo

          Df Sum of Sq    RSS     AIC
+ Gender   1    5.3699  95.35 -82.940
<none>                 100.72 -75.954
+ Age_num  1    0.0012 100.72 -73.956

Step:  AIC=-82.94
shannon ~ Expo + Gender

          Df Sum of Sq   RSS     AIC
<none>                 95.35 -82.940
+ Age_num  1  0.019726 95.33 -80.974

Start:  AIC=-33.65
shannon ~ 1

          Df Sum of Sq    RSS     AIC
+ Expo     1    6.8840 47.136 -42.962
<none>                 54.020 -33.647
+ Age_num  1    1.0369 52.983 -33.256
+ Gender   1    0.0162 54.004 -31.672

Step:  AIC=-42.96
shannon ~ Expo

          Df Sum of Sq    RSS     AIC
+ Age_num  1   1.53787 45.598 -43.715
<none>                 47.136 -42.962
+ Gender   1   0.02836 47.108 -41.012

Step:  AIC=-43.71
shannon ~ Expo + Age_num

         Df Sum of Sq    RSS     AIC
<none>                45.598 -43.715
+ Gender  1   0.04796 45.550 -41.802

Start:  AIC=1575.04
observed_otus ~ 1

          Df Sum of Sq     RSS    AIC
+ Expo     1    415511 1985849 1545.9
+ Gender   1     58732 2342629 1573.0
<none>                 2401361 1575.0
+ Age_num  1     13739 2387621 1576.1

Step:  AIC=1545.88
observed_otus ~ Expo

          Df Sum of Sq     RSS    AIC
+ Gender   1     71908 1913942 1541.8
<none>                 1985849 1545.9
+ Age_num  1         0 1985849 1547.9

Step:  AIC=1541.83
observed_otus ~ Expo + Gender

          Df Sum of Sq     RSS    AIC
<none>                 1913942 1541.8
+ Age_num  1    148.46 1913793 1543.8

Start:  AIC=822.68
observed_otus ~ 1

          Df Sum of Sq     RSS    AIC
+ Expo     1     62837 1571116 821.42
<none>                 1633952 822.68
+ Gender   1      4528 1629424 824.45
+ Age_num  1      3480 1630472 824.50

Step:  AIC=821.42
observed_otus ~ Expo

          Df Sum of Sq     RSS    AIC
<none>                 1571116 821.42
+ Age_num  1    6408.9 1564707 823.08
+ Gender   1    4018.0 1567098 823.21

Start:  AIC=912.08
faith_pd ~ 1

          Df Sum of Sq   RSS    AIC
+ Expo     1   17233.1 24922 827.88
+ Age_num  1    1860.3 40295 906.67
+ Gender   1    1403.9 40751 908.52
<none>                 42155 912.08

Step:  AIC=827.88
faith_pd ~ Expo

          Df Sum of Sq   RSS    AIC
+ Gender   1   1705.46 23217 818.25
+ Age_num  1    532.68 24389 826.33
<none>                 24922 827.88

Step:  AIC=818.25
faith_pd ~ Expo + Gender

          Df Sum of Sq   RSS    AIC
+ Age_num  1    588.06 22629 816.04
<none>                 23217 818.25

Step:  AIC=816.04
faith_pd ~ Expo + Gender + Age_num

Start:  AIC=665.58
faith_pd ~ 1

          Df Sum of Sq   RSS    AIC
+ Expo     1    7975.2 26015 636.30
+ Age_num  1     789.1 33201 664.84
+ Gender   1     641.0 33349 665.36
<none>                 33990 665.58

Step:  AIC=636.3
faith_pd ~ Expo

          Df Sum of Sq   RSS    AIC
+ Gender   1    799.88 25215 634.64
+ Age_num  1    536.00 25479 635.86
<none>                 26015 636.30

Step:  AIC=634.64
faith_pd ~ Expo + Gender

          Df Sum of Sq   RSS    AIC
+ Age_num  1    545.27 24670 634.09
<none>                 25215 634.64

Step:  AIC=634.09
faith_pd ~ Expo + Gender + Age_num

Start:  AIC=-916.45
pielou_e ~ 1

          Df Sum of Sq     RSS     AIC
+ Expo     1 0.0310501 0.57520 -923.08
+ Age_num  1 0.0218545 0.58440 -920.47
<none>                 0.60625 -916.45
+ Gender   1 0.0044785 0.60178 -915.67

Step:  AIC=-923.08
pielou_e ~ Expo

          Df Sum of Sq     RSS     AIC
+ Age_num  1 0.0148881 0.56032 -925.38
<none>                 0.57520 -923.08
+ Gender   1 0.0051937 0.57001 -922.56

Step:  AIC=-925.38
pielou_e ~ Expo + Age_num

         Df Sum of Sq     RSS     AIC
<none>                0.56032 -925.38
+ Gender  1 0.0057042 0.55461 -925.05

Start:  AIC=-538.22
pielou_e ~ 1

          Df Sum of Sq    RSS     AIC
+ Age_num  1  0.066373 1.0896 -543.14
+ Expo     1  0.020103 1.1358 -538.27
<none>                 1.1559 -538.22
+ Gender   1  0.007772 1.1482 -537.01

Step:  AIC=-543.14
pielou_e ~ Age_num

         Df Sum of Sq    RSS     AIC
<none>                1.0896 -543.14
+ Expo    1 0.0162888 1.0733 -542.90
+ Gender  1 0.0081806 1.0814 -542.02

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = regforward.out, seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 44.999  6.518

Meaningful coefficients of the linear terms:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.8078937  0.0153066  52.781  < 2e-16 ***
Age_num     -0.0027632  0.0008392  -3.293  0.00132 ** 
U1.Age_num   0.0085514  0.0041917   2.040       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0955 on 113 degrees of freedom
Multiple R-Squared: 0.1085,  Adjusted R-squared: 0.08479 

Convergence attained in 1 iter. (rel. change 3.4462e-07)

    Davies' test for a change in the slope

data:  formula = pielou_e ~ Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 46.889, n.points = 9, p-value = 0.09208
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 44.9993     32.0859    57.9127
$Age_num
             Est.    St.Err. t value  CI(95%).l  CI(95%).u
slope1 -0.0027632 0.00083918 -3.2927 -0.0044258 -0.0011006
slope2  0.0057882 0.00410680  1.4094 -0.0023482  0.0139250

Start:  AIC=-76.15
shannon ~ 1

          Df Sum of Sq     RSS      AIC
+ Expo     1   19.5685  82.264 -109.149
+ Age_num  1    4.6701  97.162  -81.851
+ Gender   1    2.0289  99.804  -77.452
<none>                 101.833  -76.152

Step:  AIC=-109.15
shannon ~ Expo

          Df Sum of Sq    RSS     AIC
+ Gender   1    2.4132 79.851 -112.03
+ Age_num  1    2.2321 80.032 -111.66
<none>                 82.264 -109.15

Step:  AIC=-112.03
shannon ~ Expo + Gender

          Df Sum of Sq    RSS     AIC
+ Age_num  1    2.3663 77.484 -114.97
<none>                 79.851 -112.03

Step:  AIC=-114.97
shannon ~ Expo + Gender + Age_num

Start:  AIC=33.32
shannon ~ 1

          Df Sum of Sq    RSS    AIC
+ Expo     1   13.5753 139.34 24.446
+ Age_num  1    8.3542 144.56 28.750
<none>                 152.92 33.323
+ Gender   1    2.1475 150.77 33.668

Step:  AIC=24.45
shannon ~ Expo

          Df Sum of Sq    RSS    AIC
+ Age_num  1    7.2321 132.11 20.210
+ Gender   1    2.5217 136.82 24.309
<none>                 139.34 24.446

Step:  AIC=20.21
shannon ~ Expo + Age_num

         Df Sum of Sq    RSS    AIC
+ Gender  1    2.5822 129.53 19.901
<none>                132.11 20.210

Step:  AIC=19.9
shannon ~ Expo + Age_num + Gender

Start:  AIC=1592.59
observed_otus ~ 1

          Df Sum of Sq     RSS    AIC
+ Expo     1    938567 1734140 1523.7
+ Age_num  1    106387 2566320 1587.9
+ Gender   1     62096 2610611 1590.7
<none>                 2672707 1592.6

Step:  AIC=1523.65
observed_otus ~ Expo

          Df Sum of Sq     RSS    AIC
+ Gender   1     76969 1657171 1518.2
+ Age_num  1     31788 1702352 1522.6
<none>                 1734140 1523.7

Step:  AIC=1518.2
observed_otus ~ Expo + Gender

          Df Sum of Sq     RSS    AIC
+ Age_num  1     34657 1622514 1516.7
<none>                 1657171 1518.2

Step:  AIC=1516.74
observed_otus ~ Expo + Gender + Age_num

Start:  AIC=1182.19
observed_otus ~ 1

          Df Sum of Sq     RSS    AIC
+ Expo     1    534368 2277124 1159.5
+ Age_num  1     71038 2740454 1181.2
+ Gender   1     57098 2754394 1181.8
<none>                 2811492 1182.2

Step:  AIC=1159.52
observed_otus ~ Expo

          Df Sum of Sq     RSS    AIC
+ Gender   1     69292 2207832 1157.9
+ Age_num  1     51145 2225979 1158.9
<none>                 2277124 1159.5

Step:  AIC=1157.91
observed_otus ~ Expo + Gender

          Df Sum of Sq     RSS    AIC
+ Age_num  1     51988 2155844 1157.1
<none>                 2207832 1157.9

Step:  AIC=1157.12
observed_otus ~ Expo + Gender + Age_num

colnames(model_tbls) = c("Body_Site", "Indice", "Year", "Model", "Model.type")
model_tbls$Model <- gsub(" + psi1.Age_num", "", model_tbls$Model, fixed = T) %>% paste0(model_tbls$Indice, " ~ ", .)

Graphs

Indices <- colnames(alpha[c(2:5)])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Nose", "Mouth", "Right_Arm", "Right_Hand")

dat_new = complete(cbind(Age_num = seq(1, 60, 1), 
                         Expo = c("Low", "Medium"), 
                         Gender = c("F", "M")) %>% 
                       as.data.frame(), Age_num, Expo, Gender) %>% 
    mutate(Age_num = as.numeric(Age_num))

dat_pred = list()
for (bb in BodySites){
    for (aa in Indices){
        for(yy in Yrs){
            dat_pred[[bb]][[aa]][[yy]] = cbind(dat_new, predict(models[[bb]][[aa]][[yy]], newdata = dat_new, interval = "confidence"))
            
        }
    }
}
for (yy in Yrs){
    for (aa in Indices){
        #aa = "faith_pd"
        p_list = list()
        for (bb in BodySites){
            #bb = "Feces"
            dat_ = dat_alpha %>% dplyr::filter(Body_Site==bb, Year==yy, SampleGroup=="Villagers", Age_num>=1) %>% mutate(Expo = as.character(Expo))
            dat_pred_ = dat_pred[[bb]][[aa]][[yy]]
            mod = model_tbls %>% dplyr::filter(Body_Site==bb, Year==yy, Indice==aa) %>% pull(Model)
            p.expo = (models[[bb]][[aa]][[yy]] %>% car::Anova(type = "3") %>% as.data.frame())["Expo", "Pr(>F)"]
            p.expo = ifelse(p.expo<0.001, "P.Expo<0.001", paste0("P.Expo=", round(p.expo, 3)))
            p = ggplot(dat_, aes_string(x = "Age_num", y = aa, color = "Expo", shape = "Gender")) +
                geom_point(size = 1.5, alpha = 0.6) +
                geom_line(data = dat_pred_, aes(y = fit, linetype = Gender), size = 1) +
                geom_ribbon(data = dat_pred_, aes(y = fit, ymin = lwr, ymax = upr, fill = Expo, alpha = Gender), color = NA, show.legend = F) +
                annotate("label", x = Inf, y = Inf, label = paste(stringi::stri_wrap(c(mod, p.expo),width = 50), collapse = "\n"), hjust = 1, vjust = 1, size = 9/.pt) +
                scale_color_manual(values = Expo_color, limits = force) +
                scale_fill_manual(values = Expo_color, limits = force) +
                scale_shape_manual(values = c(16, 1)) +
                scale_linetype_manual(values = c(1, 2)) +
                scale_alpha_manual(values = c(0.4, 0.2)) +
                theme_jw(base_size = 14) + theme(aspect.ratio = 0.7, panel.background = element_rect(color = "black"), plot.title = element_text(hjust = 0, vjust = 1), plot.title.position = "panel", legend.key.width = unit(1,"cm"), legend.margin = margin(rep(-3, 4))) +
                labs(color = "Samples points", linetype = "Regression lines", shape = "", x = "Age (Years)") +
                guides(color = guide_legend(override.aes = list(alpha = 1, size = 1.5, linetype = c(NA, NA)), order = 1), shape = guide_legend(override.aes = list(size = 1.5, alpha = 1), order = 2), linetype = guide_legend(order = 3)) +
                ggtitle(bb)
            p_list[[bb]] = p
        }
        p_f = p_list[["Feces"]] + p_list[["Mouth"]] + guide_area() + p_list[["Nose"]] + p_list[["Right_Arm"]] + p_list[["Right_Hand"]] + plot_layout(guides = "collect", nrow = 2) + plot_annotation(title = yy, tag_levels = 'A') & theme(plot.tag = element_text(size = 10, face = "bold"))
        print(p_f)
        #set_panel_size(p = p_f, g = patchworkGrob(p_f), file = paste0( "../output/figures/alpha_regression_", aa, "_", yy, ".pdf"), width =  unit(3, "in"), height = unit(2.1, "in"))
    }
}